## Austin R. Nelson, Texas Law Review
## Evaluating the Evaluator (Nov. 2020)

# Clear environment
remove(list=ls())

# Packages from library
library(readxl)
library(MASS)
library(stargazer)
library(brant)
library(ggplot2)

# Open data
ABA <- read_excel(#PUT DATA LOCATION HERE#)
View(ABA)
attach(ABA)

# Summarize data
summary(ABA)
mean(rating)
mean(prestige)
sd(rating)
sd(bush)
sd(obama)
sd(trump)
sd(male)
sd(fedjudge)
sd(appclerk)
sd(prestige)
sd(yrslawyer)

# Obama-Trump models
OTm1 <- polr(as.factor(rating) ~ obama + trump, method = "logistic", Hess = TRUE)
summary(OTm1)

OTm2 <- polr(as.factor(rating) ~ obama + trump + male + fedjudge + appclerk, method = "logistic", Hess = TRUE)
summary(OTm2)

OTm3 <- polr(as.factor(rating) ~ obama + trump + male + fedjudge + appclerk + prestige + yrslawyer, method = "logistic", Hess = TRUE)
summary(OTm3)

stargazer(OTm1, OTm2, OTm3, type = "text")

# Brant test of parallel regression assumption
# See Brant, R. 1990. "Assessing proportionality in the proportional odds model for ordinal logistic regression." Biometrics 46: 1171-78.
brant(OTm1)
brant(OTm2)
brant(OTm3)

# Calculate odds ratios and confidence intervals
exp(coef(OTm1))
exp(coef(OTm2))
exp(coef(OTm3))
exp(confint(OTm1))
exp(confint(OTm2))
exp(confint(OTm3))

cbind(exp(coef(OTm3)), exp(confint(OTm3))) # Table for model 3

# Forest plot
FPlabel <- factor(c("Obama","Trump","Male","Fed. Judge","App. Clerk","Law School","Yrs. Lawyer"), levels=c("Yrs. Lawyer","Law School","App. Clerk","Fed. Judge","Male","Trump","Obama"))
FPmean  <- c(exp(coef(OTm3))) 
FPlower <- c(exp(confint(OTm3))[,1])
FPupper <- c(exp(confint(OTm3))[,2])

FPdf <- data.frame(FPlabel, FPmean, FPlower, FPupper)
View(FPdf)

fp <- ggplot(data=FPdf, aes(x=FPlabel, y=FPmean, ymin=FPlower, ymax=FPupper)) +
  geom_pointrange() + 
  geom_hline(yintercept=1, lty=2) +  # add a dotted line at x=1 after flip
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("Odds Ratio (95% Confidence Interval)") +
  theme_bw() + # use a white background
  scale_y_continuous(trans="log",breaks=c(0.5,1,2,4,8,16)) +
  theme(text = element_text(size = 14, family = "serif")) +
  ggtitle("Figure 1. Odds Ratios for ABA Judicial Ratings, 2001-2020")
print(fp)

# MARGINGAL EFFECTS PLOT (WQ/Q/NQ)
pred.data <- data.frame(obama=mean(obama), trump=mean(trump), male=median(male), fedjudge=median(fedjudge), appclerk=median(appclerk), prestige=mean(prestige), yrslawyer=seq(min(yrslawyer), max(yrslawyer), by=1))
pred.data$id <- row.names(pred.data)
pred.data$pred0 <- predict(OTm3, pred.data, type="probs")[,1]
pred.data$pred1 <- predict(OTm3, pred.data, type="probs")[,2]
pred.data$pred2 <- predict(OTm3, pred.data, type="probs")[,3]

View(pred.data)
attach(pred.data)
plot(x = yrslawyer, y = pred2, ylim = c(0,1), pch = 19, col = "#008200", lwd = 2, main = "Figure 2. Probability of ABA Ratings", xlab = "Years as a Lawyer", ylab = "Probability", font.lab = "2", sub = "", font.sub = "3", family = "serif")
grid(col = "lightgray", lty = "solid")
points(x = yrslawyer, y = pred2, pch = 19, col = "#008200", lwd = 2)
points(x = yrslawyer, y = pred1, pch = 15, col = "black", lwd = 2)
points(x = yrslawyer, y = pred0, pch = 17, col = "red", lwd = 2)
lines(x = yrslawyer, y = pred2, col = "#008200", lwd = 1)
lines(x = yrslawyer, y = pred1, col = "black", lwd = 1)
lines(x = yrslawyer, y = pred0, col = "red", lwd = 1)
legend("right", cex = 1.0, c("Well Qualified", "Qualified", "Not Qualified"), pch = c(19,15,17), col = c("#008200","black","red"))

# MARGINGAL EFFECTS PLOT (Bush/Obama/Trump)
attach(ABA)
pred.bush <- data.frame(obama=0, trump=0, male=median(male), fedjudge=median(fedjudge), appclerk=median(appclerk), prestige=mean(prestige), yrslawyer=seq(min(yrslawyer), max(yrslawyer), by=1))
pred.bush$id <- row.names(pred.bush)
pred.bush$pred2 <- predict(OTm3, pred.bush, type="probs")[,3]
View(pred.bush)
attach(pred.bush)
plot(x = yrslawyer, y = pred2, ylim = c(0,1), pch = 19, col = "pink", lwd = 2, main = "Figure 3. Well Qualified Ratings by Nominating President", xlab = "Years as a Lawyer", ylab = "Probability", font.lab = "2", sub = "", font.sub = "3", family = "serif")
grid(col = "lightgray", lty = "solid")
points(x = yrslawyer, y = pred2, pch = 19, col = "pink", lwd = 2)
lines(x = yrslawyer, y = pred2, col = "pink", lwd = 1)

attach(ABA)
pred.obama <- data.frame(obama=1, trump=0, male=median(male), fedjudge=median(fedjudge), appclerk=median(appclerk), prestige=mean(prestige), yrslawyer=seq(min(yrslawyer), max(yrslawyer), by=1))
pred.obama$id <- row.names(pred.obama)
pred.obama$pred2 <- predict(OTm3, pred.obama, type="probs")[,3]
View(pred.obama)
attach(pred.obama)
points(x = yrslawyer, y = pred2, pch = 19, col = "blue", lwd = 1)
lines(x = yrslawyer, y = pred2, col = "blue", lwd = 1)

attach(ABA)
pred.trump <- data.frame(obama=0, trump=1, male=median(male), fedjudge=median(fedjudge), appclerk=median(appclerk), prestige=mean(prestige), yrslawyer=seq(min(yrslawyer), max(yrslawyer), by=1))
pred.trump$id <- row.names(pred.trump)
pred.trump$pred2 <- predict(OTm3, pred.trump, type="probs")[,3]
View(pred.trump)
attach(pred.trump)
points(x = yrslawyer, y = pred2, pch = 19, col = "red", lwd = 1)
lines(x = yrslawyer, y = pred2, col = "red", lwd = 1)

legend("bottomright", cex = 1.0, c("Obama", "Trump", "Bush"), fill=c("blue", "red", "pink"))
